% This code produces the figures with the response of markups
% by Raphael Schoenle, February 18, 2016
% updated September 6, 2016

% six cases with parameterizations:
% one basecase and one parameter variation (phi_pi=2.5)
clc; 
clear all;

cd /Users/schoenle/Dropbox/Networks/code/gensys/outputs_MP_paper_main_new

%number of sectors
n_sectors = 350;
n_sectors1 = 58;

%number of periods in IRF
n_periods = 30;
n_periods_idio = 30;

%number of cases
num_case = 18;

%pre-definition of variables
% third dimension: number of shocks: monetary, technology, idiosyncratic
c_imp = zeros(n_periods, num_case, 3);
inf_imp = zeros(n_periods, num_case, 3);
% fourth dimension: number of sectors, max set at 350
mc_imp = zeros(n_periods, num_case, 3, n_sectors);
markup_imp = zeros(n_periods, num_case, 3, n_sectors);
pos_akshock = zeros(n_sectors,1);
pos_akshock1 = zeros(n_sectors1,1);

%Parameters values
params.rho_a = 0.9;         % persistence agg tech shocks
params.rho_mu = 0.9;         % persistence monetary shock
params.rho_k = 0;         % persistence idiosync tech shock
params.beta = 0.9975;      % discount factor
params.theta = 6;          % elast of subst within sectors

% parameter values to loop over
% a. eta = 2, phi = 2, sigma = 1, phi_pi = 1.24, phi_c = 0.33/12 (base case
% in paper)
% b. eta = 2, phi = 2, sigma = 1, phi_pi = 2.5, phi_c = 0.33/12 same as a.
% but with higher response to inflation in Taylor rule

% II. delta = 0.5
% a. - g. as above

% parameter sets
delta_set = [0.5 0.5];
eta_set = [2 2];
phi_set = [2 2];
sigma_set = [1 1];
phi_pi_set = [1.24 2.5];
phi_c_set = [0.33/12 0.33/12];
lambda1_set =(1+phi_pi_set)./(sigma_set+ phi_c_set);
lambda2_set = 1./(sigma_set + phi_c_set); 
lambda3_set = sigma_set./(sigma_set + phi_c_set);
lambda4_set = lambda2_set;
lambda5_set = phi_pi_set./(sigma_set + phi_c_set);
clear l1 l2 l3 l4 l5;
n_parameters = length(delta_set);

% Position indices for shocks
pos_mushock = 1;
pos_ashock = 2;
for i = 1:n_sectors
    pos_akshock(i) = 2 + i;
end
for i = 1:n_sectors1
    pos_akshock1(i) = 2 + i;
end

% sectoral weights in consumption
wc_het = xlsread('/Users/schoenle/Dropbox/Networks/IO_tables/gdpshares_2002_valad_concorded.xlsx');
Omega_C_het = wc_het(:,2)./sum(wc_het(:,2));
Omega_C_hom = 1/n_sectors*ones(size(Omega_C_het));

wc_58_het = xlsread('/Users/schoenle/Dropbox/Networks/IO_tables/gdpshares_58_2002_valad_concorded_new.xlsx');
Omega_C_58_het = wc_58_het(:,2);
Omega_C_58_hom = 1/n_sectors1*ones(size(Omega_C_58_het));

% sectoral Calvo
% 350 sectors first
zzz = xlsread('/Users/schoenle/Dropbox/Networks/IO_tables/secfreq_new.xlsx');
FPA_het  = zzz(:,2);
FPA_hom = mean(FPA_het)*ones(size(FPA_het));
FPA_flex = 0.99*ones(size(FPA_het));
FPA_sticky = 0.01*ones(size(FPA_het));
% Omega_C weighted mean
FPA_hom_weighted = Omega_C_het'*FPA_het*ones(size(FPA_het));

% 58 sectors
FPA_58_het = xlsread('/Users/schoenle/Dropbox/Networks/IO_tables/secfreq_58.xlsx');
FPA_58_het = (mean(FPA_het)-mean(FPA_58_het))+FPA_58_het;
FPA_58_hom = mean(FPA_het)*ones(size(FPA_58_het));
FPA_58_flex =  0.99*ones(size(FPA_58_het));
FPA_58_sticky =  0.01*ones(size(FPA_58_het));
FPA_58_hom_weighted = Omega_C_58_het'*FPA_58_het*ones(size(FPA_58_het));

%sectoral weights 
Omega_het = xlsread('/Users/schoenle/Dropbox/Networks/IO_tables/IO_links.xlsx')';
Omega_het = (Omega_het./repmat(sum(Omega_het,2),1,n_sectors));
Omega_hom = 1/n_sectors*ones(size(Omega_het)); 
Omega_C_as_IO = repmat(Omega_C_het',n_sectors,1);

Omega_58_het = xlsread('/Users/schoenle/Dropbox/Networks/IO_tables/IO_links_58.xlsx')';
Omega_58_het = Omega_58_het./repmat(sum(Omega_58_het,2),1,n_sectors1);
Omega_58_hom = 1/n_sectors1*ones(size(Omega_58_het)); 
Omega_58_C_as_IO = repmat(Omega_C_58_het',n_sectors1,1);

%% 350 Sectors
%Case1: flex-price Calvo, homogeneous size, homogeneous I/O
FPA(:,:,1) = FPA_flex;
Omega_C(:,:,1) = Omega_C_hom;
Omega(:,:,1) = Omega_hom;

%Case2: homogeneous weighted mean Calvo, homogeneous size, homogeneous I/O
FPA(:,:,2) = FPA_hom_weighted;
Omega_C(:,:,2) = Omega_C_hom;
Omega(:,:,2) = Omega_hom;

%Case3: heterogeneous Calvo, homogeneous size, homogeneous I/O
FPA(:,:,3) = FPA_het;
Omega_C(:,:,3) = Omega_C_hom;
Omega(:,:,3) = Omega_hom;

%Case4: heterogeneous Calvo, heterogeneous size, heterogeneous
%I/O based on size weights
FPA(:,:,4) = FPA_het;
Omega_C(:,:,4) = Omega_C_het;
Omega(:,:,4) = Omega_C_as_IO;

%Case5: heterogeneous Calvo, heterogeneous size, homogeneous I/O
FPA(:,:,5) = FPA_het;
Omega_C(:,:,5) = Omega_C_het;
Omega(:,:,5) = Omega_hom;

%Case6: heterogeneous Calvo, heterogeneous size, heterogeneous I/O
FPA(:,:,6) = FPA_het;
Omega_C(:,:,6) = Omega_C_het;
Omega(:,:,6) = Omega_het;


%% 58 Sectors
%Case1: flex-price Calvo, homogeneous size, homogeneous I/O
FPA_58(:,:,1) = FPA_58_flex;
Omega_C_58(:,:,1) = Omega_C_58_hom;
Omega_58(:,:,1) = Omega_58_hom;

%Case2: homogeneous weighted mean Calvo, homogeneous size, homogeneous I/O
FPA_58(:,:,2) = FPA_58_hom_weighted;
Omega_C_58(:,:,2) = Omega_C_58_hom;
Omega_58(:,:,2) = Omega_58_hom;

%Case3: heterogeneous Calvo, homogeneous size, homogeneous I/O
FPA_58(:,:,3) = FPA_58_het;
Omega_C_58(:,:,3) = Omega_C_58_hom;
Omega_58(:,:,3) = Omega_58_hom;

%Case4: heterogeneous Calvo, heterogeneous size, heterogeneous
%I/O based on size weights
FPA_58(:,:,4) = FPA_58_het;
Omega_C_58(:,:,4) = Omega_C_58_het;
Omega_58(:,:,4) = Omega_58_C_as_IO;

%Case5: heterogeneous Calvo, heterogeneous size, homogeneous I/O
FPA_58(:,:,5) = FPA_58_het;
Omega_C_58(:,:,5) = Omega_C_58_het;
Omega_58(:,:,5) = Omega_58_hom;

%Case6: heterogeneous Calvo, heterogeneous size, heterogeneous I/O
FPA_58(:,:,6) = FPA_58_het;
Omega_C_58(:,:,6) = Omega_C_58_het;
Omega_58(:,:,6) = Omega_58_het;

% index for the most/least flexible sector
for j=1:num_case/3
	i_max_FPA_58(j) = find(FPA_58(:,:,j)==max(FPA_58(:,:,j)),1);
	i_min_FPA_58(j)= find(FPA_58(:,:,j)==min(FPA_58(:,:,j)),1);
	i_max_FPA_350(j) = find(FPA(:,:,j)==max(FPA(:,:,j)),1);
	i_min_FPA_350(j) = find(FPA(:,:,j)==min(FPA(:,:,j)),1);
end

% index for the largest/smallest sector
for j=1:num_case/3
	i_max_Omega_C_58(j) = find(Omega_C_58(:,:,j)==max(Omega_C_58(:,:,j)),1);
	i_min_Omega_C_58(j)= find(Omega_C_58(:,:,j)==min(Omega_C_58(:,:,j)),1);
	i_max_Omega_C_350(j) = find(Omega_C(:,:,j)==max(Omega_C(:,:,j)),1);
	i_min_Omega_C_350(j) = find(Omega_C(:,:,j)==min(Omega_C(:,:,j)),1);
end


% outer loop over the parameter values, index k
for k = 1:1
    params.delta = delta_set(k);
    params.eta = eta_set(k);
    params.frisch = phi_set(k);
    params.sigma = sigma_set(k);
    params.phi_pi = phi_pi_set(k);
    params.phi_c = phi_c_set(k);
    params.lambda1 = lambda1_set(k);
    params.lambda2 = lambda2_set(k);
    params.lambda3 = lambda3_set(k);
    params.lambda4 = lambda4_set(k);
    params.lambda5 = lambda5_set(k);

    % inner loop over the cases, index i
    % 1-6: 350 sectors
    % 7-12: 58 sectors
    
    for i = 1:12
    i
    tic
    if i < num_case/3+1
        % 350-sector cases	
        %%Generate Gensys matrices
        [G0, G1, C, PSI, PI, pos_c, pos_pc, pos_mc, pos_yk, pos_p, pos_ps]  = Gensys_matrix(n_sectors, FPA(:,:,i), Omega_C(:,:,i), Omega(:,:,i), params);
    toc
        %%Call Gensys
        [G1, C, imp, fmat, fwt, ywt, gev, eu, loose] = gensys(G0, G1, C, PSI, PI);
    toc
        %%Generate IRFs
        %%%Monetary policy shocks
        % now include MC as output and markups (P_k/MC_k)
        yt_imp_mu = Gensys_IRF(G1, n_sectors, n_periods, pos_mushock, imp);
        output = [1:n_periods; yt_imp_mu(pos_c,2:end); yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1); yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end) ),n_sectors,1);   yt_imp_mu(pos_p,2:end) -   yt_imp_mu(pos_mc,2:end)   ]';
        autocorr_c = autocorr(output(:,end-1),1);
        autocorr_inf = autocorr(output(:,end),1);
        autocorr_mc = autocorr(mean([yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end)),n_sectors,1)]',2),1);        
        % output matrix: initial C, Pi, mean(MC) cumulative C, Pi, mean(MC)...
        % AR(1) of C, Pi, MC + IRFs (C, Pi, MC)
        output_matrix = [output(:,1) repmat([output(1,2:3)...
            mean([yt_imp_mu(pos_mc,2)-repmat((yt_imp_mu(pos_pc,2)),n_sectors,1)]',2)...
        sum(output(:,2:3)) sum(mean([yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end)),n_sectors,1)]',2))...    
        autocorr_c(2) autocorr_inf(2) autocorr_mc(2)],n_periods,1) output(:,2:end)];   
        
        % cumulative responses - sort and get sector indices
        [output_top_bottom, I] = sort([sum(yt_imp_mu(pos_yk,2:end),2)]);
        top_bottom = [I(1:10) , I(end-9:end)];
        
        % save output
        filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_350.csv'];
        fid = fopen(filename, 'w');
        fprintf(fid,'Period\tConsumption_impact_irf_mu\tInflation_impact_irf_mu\tMC_mean_impact_irf_mu\tConsumption_cum_irf_mu\tInflation_cum_irf_mu\tMC_mean_cum_irf_mu\tConsumption_persistence_mu\tInflation_persistence_mu\tMC_persistence_mu\tConsumption_mu\tInflation_mu\tMC_mu \t Pk/MCk \n');
        fclose(fid);
        dlmwrite(filename, (output_matrix), '-append', 'precision', '%.10f', 'delimiter', '\t');

        mc_imp(:,i,1,:) = [yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end)),n_sectors,1)]';
        markup_imp(:,i,1,:) = [ yt_imp_mu(pos_p,2:end) -   yt_imp_mu(pos_mc,2:end)  ]';
        
    elseif i>num_case/3 && i<=num_case/3*2
        clear yt_imp_a_idio
        [G0, G1, C, PSI, PI, pos_c, pos_pc, pos_mc, pos_yk, pos_p, pos_ps]  = Gensys_matrix(n_sectors1, FPA_58(:,:,i-num_case/3), Omega_C_58(:,:,i-num_case/3), Omega_58(:,:,i-num_case/3), params);    
    toc
        %%Call Gensys
        [G1, C, imp, fmat, fwt, ywt, gev, eu, loose] = gensys(G0, G1, C, PSI, PI);
    toc
        %%Generate IRFs
        % now include MC as output
        yt_imp_mu = Gensys_IRF(G1, n_sectors1, n_periods, pos_mushock, imp);
        output = [1:n_periods; yt_imp_mu(pos_c,2:end); yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1); yt_imp_mu(pos_mc,2:end)-repmat(yt_imp_mu(pos_pc,2:end),n_sectors1,1);  yt_imp_mu(pos_p,2:end) -   yt_imp_mu(pos_mc,2:end)  ]';
        autocorr_c = autocorr(output(:,2),1);
        autocorr_inf = autocorr(output(:,3),1);
        autocorr_mc = autocorr(mean([yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end)),n_sectors1,1)]',2),1);        
        % output matrix: initial C, Pi, mean(MC) cumulative C, Pi, mean(MC)...
        % AR(1) of C, Pi, MC + IRFs (C, Pi, MC)
        output_matrix = [output(:,1) repmat([output(1,2:3)...
            mean([yt_imp_mu(pos_mc,2)-repmat((yt_imp_mu(pos_pc,2)),n_sectors1,1)]',2)...
        sum(output(:,2:3)) sum(mean([yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end)),n_sectors1,1)]',2))...    
        autocorr_c(2) autocorr_inf(2) autocorr_mc(2)],n_periods,1) output(:,2:end)];   
        % cumulative responses - sort and get sector indices
        [output_top_bottom, I] = sort([sum(yt_imp_mu(pos_yk,2:end),2)]);
        top_bottom = [I(1:10) , I(end-9:end)];
        
        % save output
        filename = ['Gensys_mu_output_main_case_' num2str(i-num_case/3) 'params' num2str(k) '_58.csv'];
        fid = fopen(filename, 'w');
        fprintf(fid,'Period\tConsumption_impact_irf_mu\tInflation_impact_irf_mu\tMC_mean_impact_irf_mu\tConsumption_cum_irf_mu\tInflation_cum_irf_mu\tMC_mean_cum_irf_mu\tConsumption_persistence_mu\tInflation_persistence_mu\tMC_persistence_mu\tConsumption_mu\tInflation_mu\tMC_mu \t Pk/MCk\n');
        fclose(fid);
        dlmwrite(filename, (output_matrix), '-append', 'precision', '%.10f', 'delimiter', '\t');

        mc_imp(:,i,1,1:n_sectors1) = [yt_imp_mu(pos_mc,2:end)-repmat((yt_imp_mu(pos_pc,2:end)),n_sectors1,1)]';

        markup_imp(:,i,1,1:n_sectors1) = [ yt_imp_mu(pos_p,2:end) -   yt_imp_mu(pos_mc,2:end)  ]';

    end

    %%Save Case output in output matrix
    c_imp(:,i,1) = yt_imp_mu(pos_c,2:end)';
    inf_imp(:,i,1) = (yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1))';
    
    toc
    end
    


cd /Users/schoenle/Dropbox/Networks/code/gensys/outputs_MP_paper_main_new


    
    
    
    
    
    
    
    
    % case by case plot of markups:
    % mean, median, and percentiles: 1, 5, 25, 75, 95, 99
    figure
    
    %('units','normalized','position',[.1 0 .4 1])
    i = 2;
    % 350 goods, mean markup
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),1,4)','-b','linewidth',2)
    hold on
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),5,4)','-bx','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),25,4),'--b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),50,4),'-.b','linewidth',4)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),75,4),'-.b+','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),95,4),'-.b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),99,4)',':b','linewidth',2)
    %title(['Markups, Case ', num2str(i-num_case/3)],'Interpreter','latex')
    title('Panel A')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    legend1=legend('1st Percentile','5th Percentile','25th Percentile','Median','75th Percentile','95th Percentile','99th Percentile','best')
    axis('square')

    fig58 = gcf
    fig58.PaperUnits = 'normalized';
    fig58.PaperPosition = [.1 0 .8 1];
    fig58.PaperPositionMode = 'manual';
    filename = ['Combined_IRFs_markups_percentiles_sectors_panelA_350.pdf'];
    saveas(fig58,filename);
    
    


    % case by case plot of markups:
    % mean, median, and percentiles: 1, 5, 25, 75, 95, 99
    figure
    
    %('units','normalized','position',[.1 0 .4 1])
    i = num_case/3+2;
    % 58 goods, mean markup
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),1,4)','-b','linewidth',2)
    hold on
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),5,4)','-bx','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),25,4),'--b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),50,4),'-.b','linewidth',4)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),75,4),'-.b+','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),95,4),'-.b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),99,4)',':b','linewidth',2)
    %title(['Markups, Case ', num2str(i-num_case/3)],'Interpreter','latex')
    title('Panel A')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    legend1=legend('1st Percentile','5th Percentile','25th Percentile','Median','75th Percentile','95th Percentile','99th Percentile','best')
    axis('square')

    fig58 = gcf
    fig58.PaperUnits = 'normalized';
    fig58.PaperPosition = [.1 0 .8 1];
    fig58.PaperPositionMode = 'manual';
    filename = ['Combined_IRFs_markups_percentiles_sectors_panelA_58.pdf'];
    saveas(fig58,filename);
    
    figure
    i = num_case/3+6;
    % 58 goods, mean markup
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),1,4)','-b','linewidth',2)
    hold on
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),5,4)','-bx','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),25,4),'--b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),50,4),'-.b','linewidth',4)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),75,4),'-.b+','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),95,4),'-.b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors1),99,4)',':b','linewidth',2)
    %title(['Markups, Case ', num2str(i-num_case/3)],'Interpreter','latex')
    title('Panel B')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
%    legend1=legend('1st Percentile','5th Percentile','25th Percentile','Median','75th Percentile','95th Percentile','99th Percentile','best')
    axis('square')

    fig58 = gcf
    fig58.PaperUnits = 'normalized';
    fig58.PaperPosition = [.1 0 .8 1];
    fig58.PaperPositionMode = 'manual';
    filename = ['Combined_IRFs_markups_percentiles_sectors_panelB_58.pdf'];
    saveas(fig58,filename);


    
    figure
    i = 6;
    % 58 goods, mean markup
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),1,4)','-b','linewidth',2)
    hold on
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),5,4)','-bx','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),25,4),'--b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),50,4),'-.b','linewidth',4)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),75,4),'-.b+','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),95,4),'-.b','linewidth',2)
    plot(1:n_periods,prctile(markup_imp(:,i,1,1:n_sectors),99,4)',':b','linewidth',2)
    %title(['Markups, Case ', num2str(i-num_case/3)],'Interpreter','latex')
    title('Panel B')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
%    legend1=legend('1st Percentile','5th Percentile','25th Percentile','Median','75th Percentile','95th Percentile','99th Percentile','best')
    axis('square')

    fig58 = gcf
    fig58.PaperUnits = 'normalized';
    fig58.PaperPosition = [.1 0 .8 1];
    fig58.PaperPositionMode = 'manual';
    filename = ['Combined_IRFs_markups_percentiles_sectors_panelB_350.pdf'];
    saveas(fig58,filename);
    
end

save('alloutput_markup.mat')